• GAMES103 - Intro to Physics-Based Animation

      (Based on Unity, C#)

      Huamin Wang (games103@style3D.com) Video; Lecture site

       

      Lecture 1 Introduction

      Graphics

      Geometry

      Rendering

      Material Scan

      Animations

      Physics-Based Animation Topics

      Hybrid Methods and Coupling - problems

       

      Lecture 2 Math Background

      Vector, Matrix and Tensor Calculus

      Vectors

      Definitions

      A geometric entity endowed with magnitude and direction

      p=[pxpypz]R3 ;o=[000]

      image-20211109004814489 image-20211109004854038

      Can also be stacked up to form a high-dim vector -> describe the state of an obj (not a geometric vector but a stacked vector)

      Arithematic

      Addition and Substraction

      (commutative)

      p±q=[px±qxpy±qypz±qz] ;p+q=q+p

      image-20211109005357876

      Example: p(t)=p+tv to represent the movement of a particle. Segment: 0<t<1 ; Ray: 0<t; Line: tR (t is an interpolant)

      image-20211109005650429image-20211109005758487

      Vector Norm

      Magnitude of a vector (length)

      Usage:

      Dot Product

      (inner product)

      <p,q>≡pq=pxqx+pyqy+pzqz=pTq=||p|| ||q||cosθ

      image-20211109104509904

      Properties:

      Example: Particle-Line Projection

      By def: s=||qo|| cosθ => s=osv¯

      s=||qo|| cosθ=||qo|| ||v||cosθ /||v||=(qo)Tv/||v||=(qo)Tv¯(normalized)

      image-20211109104953250

      Example: Plane Representation

      Check the relationship between point p and plane (s - the signed distance to the plane -> collision check / …; n - normal)

      s=(pc)Tn {>0Above the plane=0On the plane<0Below the plane

      image-20211109105816092

      Example: Particle-Sphere Collision

      ||p(t)c||2=r2(pc+tv)(pc+tv)=𝑟2(vv)t2+2(pc)vt+(pc)(pc)r2=0

      image-20211109105955443

      t is the root -> No root (no collision) / One root (tangentially) / Two roots (the first point, smaller t)

      Cross Product

      r=p×q=[pyqzpzqypzqxpxqzpxqypyqx]

      image-20211109110925156

      Properties

      Example: Trinagle Normal and Area

      image-20211109112748339

      Example: Trianle Inside / Outside Test (Co-plane)

      image-20211109113822850

      Example: Barycentric Coordinates

      12(x0p)×(x1p)n={12(x0p)×(x1p)inside12(x0p)×(x1p)outside

      Signed Areas:

      A2=12(x0p)×(x1p)nA0=12(x1p)×(x2p)nA1=12(x2p)×(x0p)nA=A0+A1+A2

      Barycentric weights of p: b0=A0/A, b1=A1/A, b2=A2/A (b0+b1+b2=1)

      Barycentric Interpolation: p=b0x0+b1x1+b2x2

      image-20211109114230534

      => Gourand Shading: Using barycentric interpolation (no longer popular)

      image-20211109114950142

      Example: Tetrahedral Volume

      image-20211109115747395

      Example: Barycentric Weight in Tetrahedra

      p splits the tetrahedron into 4 tetrahedra: V0=Vol(x3,x2,x1,p) , ….

      p inside: if and only if V0,V1,V2,V3>0

      Barycentric weights: bn=Vn/V, p=b0x0+b1x1+b2x2+b3x3

      Example: Particle-Triangle Intersection

      Matrices

      Definitions

      A=[a00a01a02a10a11a12a20a21a22]=[a0a1a2]R3×3

      Orthogonality

      An orthogonal matrix is a matrix made of orthogonal unit vectors.

      A=[a0a1a2] , such that  aiTaj={1,if i=j0,if ij
      ATA=[a0Ta1Ta2T][a0a1a2]=[a0Ta0a0Ta1a0Ta2a1Ta0a1Ta1a1Ta2a2Ta0a2Ta1a2Ta2]=I ;AT=A1

      Matrix Transformation

      Linear Solver

      Ax=b (A - Square matrix; x - Unknown to be found; b - Boundary vector)

      It’s expensive to compute A1, especially if A is large and sparse (Cannot use x=A1b)

      Direct Linear Solver

      LU factorization (Alt: Choleskey, LDLT, etc.)

      A=LU=[l00l10l11]L[un1,n1un1,nun,n]U

      First solve: (up -> down)

      Ly=b[l00l10l11][y0y1]=[b0b1]{y0=b0/l00y1=(b1l10y0)/l11...

      Then solve: (down -> up)

      Ux=y[un1,n1un1,nun,n][xn1xn]=[yn1yn]{xn=yn/un,nxn1=(yn1un1,nxn)/un1,n1

      Properties:

      Iterative Linear Solver

      x[k+1]=x[k]+αM1(bAx[k])

      (α - relaxation; M - Iterative matrix; bAx[k] - Residual error (for perfect solution -> = 0))

      Converge property: (bAx[0]=const at first)

      bAx[k+1]=bAx[k]αAM1(bAx[k])=(IαAM1)(bAx[k])=(IαAM1)k+1(bAx[0])

      So: (ρ(IαAM1) - Spectral radius (largest absolute value of eigenvalues))

      bAx[k+1]0, if ρ(IαAM1)<1

      For M, must be easy to solve, e.g. M=diag(A) (Jacobi), or M=lower(A) (Gauss-Seidel)

      Accelerate converge methods: Chebyshev, Conjugate Gradient, …

      Properties:

      Tensor Calculus

      Basic Concepts

       

      Lecture 3 Rigid Body Dynamics

      (Single rigid body: dynamics / rotation / …)

      Rigid bodies: Assume no deformations

      The goal of simulation is to update the state var. s[k] ove

      image-20211117163642294

      Translation Motion

      image-20211117215600450

      For translation motion, the state variable contains the position x and the velocity v (M - Mass; Force can be the function of pos, vel, t, …) -> Solve integral

      {v(t[1])=v(t[0])+M1t[0]t[1]f(x(t),v(t),t) dtx(t[1])=x(t[0])+t[0]t[1]v(t) dt

      Integration Methods Explained

      By def, the integral of x(t)=v(t) dt is the area.

      Type of Forces

      Translation Only Simulation

      image-20211118002125525

      Steps: (Mass M and Timestep Δt are user spec var)

      Rotation Motion

      Rotation Representations

      Rotation Represented by Matrix

      R=[r00r01r02r10r11r12r20r21r22]

      Suitable in graphics, rotation for vertices

      Not suitable for dynamics:

      Rotation Represented by Euler Angles

      Use 3 axial rotations to represent one general rotation. Each axial rotation uses an angle.

      Used in Unity. (the order is rotaion-by-Z / X / Y) Intuitive.

      Not suitable for dynamics:

      Rotation Represented by Quaternion

      Complex multiplications: In the complex system, two numbers represent a 2D point. => Quaternion: i, j, k are imaginary numbers (3D space) Four numbers represent a 3D point (with multiplication and division).

      image-20211118004837540

      Arithematic

      Let q=[s  v] be a quaternion made of 2 parts: a scalar s and a 3D vector v for ijk

      In Unity: provide multiplication, but no addition/subtraction/…; Use w, x, y, z -> s, v

      Representation

      Rotate around v by angle θ

      {q=[cosθ2v]q=1{q=[cosθ2v]v2=sin2θ2

      image-20211118011259709

      Convertible to the matrix:

      R=[s2+x2y2z22(xysz)2(xz+sy)2(xy+sz)s2x2+y2z22(yzsx)2(xzsy)2(yz+sx)s2x2y2+z2]

      Rotation Representations in Unity

      image-20211118011526881

      Rotation Motion

      image-20211118011710206 image-20211118012010411

      Use a 3D vector ω to denote angular velocity:

      The dir of ω -> the axis; The magnitude of ω -> the speed (sim to the representation of quaternion)

      Torque and Inertia

      Torque

      (Original state: ri -> Rotated: Rri, fi is a force)

      The rotational equiv of a force, describing the rotational tendency caused by a force.

      τi: perpendicular to both Rri and fi; proportional to ||Rri|| and ||fi||; proportional to sinθ (the angle of the two vectors)

      τiRri×fi

      image-20211123174242025

      Inertia

      Describes the resistance to rotational tendency caused by torque (not const)

      Left side (heavier point far away from the torque) has higher resistance (inertia) to the rotational tendency, slower rotation

      image-20211123174709048

      Ref state inertia, change with rotation (dep on pose). But no need to re-compute every time

      Iref=mi(riTri1ririT)

      (1 - 3x3 identity matrix)

      I=mi(riTRTRri1RririTRT)=mi(RriTri1RTRririTRT)=miR(riTri1ririT)RT=RIrefRT
      Use torque to represent

      image-20211118103522555

      Rigid Body Simulation

      Translational and Rotational Motion

      Rigid Body Simulation Process

      image-20211118104825451

      In Unity: No 3x3 matrices, only 4x4 (use 4x4 and set the last col / line ); Provide .inverse to inverse; …

      Implementation

      In practice, we update the same state var s={v,x,ω,q}

      Issues

       

      Lecture 4 Rigid Body Contacts

      Particle Collision Detection and Response

      Distance Functions

      Signed Distance Function

      Use a signed distance func ϕ(x) to define the distance indicating which side as well (corresponding to 0 surface)

      Examples

      image-20211123201917848

      ϕ(x)=(xp)nϕ(x)=||xc||rϕ(x)=||xp||2((xp)n)2r

      Intersection of Signed Distance Functions

      => Bool operations

      image-20211123205109282

      if ϕ0(x)<0 and ϕ1(x)<0 and ϕ2(x)<0 then inside ϕ(x)=max(ϕ0x,ϕ1(x),ϕ2(x)) // all val are negative, the max one is the closest else ϕ(x)=? // not relavent (no collision)

      Union of Signed Distance Functions

      image-20211123211143939

       

      if ϕ0(x)<0 or ϕ1(x)<0 then inside ϕ(x)min(ϕ0(x),ϕ1(x)) // approximate -> correct near outer boundary else outside ϕ(x)=min(ϕ0(x),ϕ1(x))

      -> We can consider collision detection with the union of two objects as collision detection with two separate objects

      Penalty methods

      (Implicit integration is better)

      Quadratic Penalty Method

      Check if collide - Yes -> Apply a force at the point (in the next update)

      For quadratic penalty (strength) potential, the force is linear

      image-20211123212648589

      image-20211123212726821

      Problem: Already inside -> cause artifacts

      => Add buffer help less the penetration issue (cannot strictly prevent)

      image-20211123212952453

      Problem: k too low -> buffer not work well; k too high -> too much force generated (overshooting)

      Log-Barrier Penalty Method

      Ensures that the force can be large enough, but assumes ϕ(x)<0 never happens => By adjusting Δt

      Always apply the penalty force: fρ1ϕ(x)N (ρ - Barrier strength)

      image-20211123213631306

      Problems: cannot prevent overshooting when very close to the surface; when penatration occurs, the penatration will be higher and higher (-> smaller step size -> higher costs)

      => Log-Barrier limited within a buffer as well to solve

      Frictional contacts are difficult to handle

      Impulse method

      Update the vel and pos as the collision occurs

      Changing the position:

      image-20211123223036097

      Changing the velocity: (μN - bounce coefficient, [0,1], a - frictional decay of vel)

      image-20211123224436596image-20211123224520049

      a should be minimized but not violating Coulomb’s law

      vTnew vTμTvNnew vN(1a)vTμT(1+μN)vN

      Therefore:

      amax(1μT(1+μN)vN/vTdynamic friction,0static friction)

      Can precisely control the friction effects

      Rigid Collision Detection and Response by Impulse

      Rigid Body Collision Detection

      When the body is made of many vertices, test each vertex: xix+Rri (from the mass center to the vertices)

      => detection: transverse every point if ϕ(x)<0

      Rigid Body Collision Response by Impulse

      Vertex i: (v - linear vel; ω - angular vel)

      {xix+Rri(Position)viv+ω×Rri(Velocity)

      image-20211123230332712

      But cannot modify xi and vi directly since they are not state var. (in Lec.3: 4 var are vel of mass center / pos of mass center / rotational pose / angular vel)

      Applying an impulse j at vertex i: (Δv=jM (Newton’s Law); Rri×j - torque induced by j)

      {vnewv+1Mjωnewω+I1(Rri×j)
      vinew=vnew+ωnew×Rri=v+1Mj+(ω+I1(Rri×j))×Rri=vi+1Mj(Rri)×(I1(Rri×j))

      Cross Product as a Matrix Product

      Convert the cross prod r× into a matrix prod r

      r×q=[ryqzrzqyrzqxrxqzrxqyryqx]=[0rzryrz0rxryrx0][qxqyqz]=rq

      In our case:

      vinew=vinew =vi+1Mj(Rri)I1(Rri)j

      Therefore, replace with some matrix K. Finally j can be computed with the following equations.

      vinew vi=KjK1M1(Rri)I1(Rri)

      Implementation

      image-20211124001811694

      Other details:

      Shape Matching

      Basic Idea

      Allow each vertex to have its own velocity, so it can move by itself

      image-20211124002525946

      Mathematical Formulation

      image-20211124002656804

      Want to find c and R (c - mass center): want the final rigid body (a square) close enough to the trapozoid

      (A - a matrix, not only corresponding to rotation, Ari=0 since the mass center was set to 0; E - The objective, =12c+Ariyi2)

      {c,R}=argmini12c+Rriyi2{c,A}=argmini12c+Ariyi2

      For mass center c and matrix A (Find derivatives):

      Ec=ic+Ariyi=icyi=0c=1Niyi(average)
      EA=i(c+Ariyi)riT=0A=(i(yic)riT)(iririT)1=Polar DecompositionRrotationSdeformation

      Remind: Singular value decomposition: A=UDVT (rotaion, scaling and rotation)

      image-20211109165409060

      Rotate the object back before the final rotation: A=(UVT)(VDVT)=RS (Local deformation: VDVT=S)

      image-20211124004551525

      Implementation

      Physical quantities are attached to each vertex, not to the entire body.

      image-20211124005104997

      (The function of Polar(A) is provided, ultilizing the polar deposition technique)

      Properties:

       

      Lecture 5 Physics-Based Cloth Simulation

      A Mass-Spring System

      Spring Systems

      An Ideal Spring

      Satisfies Hooke’s Law: The spring force tries to restore the rest length (k - Spring Stiffness)

      image-20211130163240282

      E(x)=12k(xL)2 ;f(x)=dEdx=k(xL)

      image-20211130163615639

      E(x)=12k(xixjL)fi=fj{fi(x)=iE=k(xixjL)xixjxixjfj(x)=jE=k(xjxiL)xjxixjij

      Multiple Springs

      The energies and forces can be simply summed up

      image-20211130164434830

      E=j=03Ej=j=03(12k(xixjL)2)fi=iE=j=03(k(xixjL)xixjxixj)

      Structures in Simulations

      Structured Spring Networks

      image-20211130164711386

      Unstructured Spring Networks

      Unstructured triangle mesh

      -> the edges into spring networks (usually in cloth simulations)

      Blue lines for bending resistance (every neighboring triangle pair)

      Triangle Mesh Representation

      Two arrays: Vertex & Triangle lines

      image-20211130165206242

      Not only edges but also an inner one (each triangle has 3 edges but there are repeated ones)

      Topological Construction

      Sort triangle edge triples: edge vertex index 0 & 1 and triangle index (index 0 < index)

      Integrators

      Explicit Integration

      Scheme

      (Notations: mi - Mass of vertex i; E - Edge list; L - Edge length list (pre-computed))

      Problem

      Overshooting: when k and/or Δt is too large

      Naive solution: Reduce Δt => Slow down the simulation

      Implicit Integration

      General Scheme (Euler Method)

      Integrate both x and v implicitly (MR3N×3N - Mass matrix (usually diagonal))

      {v[1]=v[0]+ΔtM1f[1]x[1]=x[0]+Δtv[1]or{x[1]=x[0]+Δtv[0]+Δt2M1f[1]v[1]=(x[1]x[0])/Δt

      v[1] & x[1] are unknown for the current time step => Find the x and v:

      Assume f dep only on x (homonomic): Solve the eqn (Problem: f may not be a linear function)

      x[1]=x[0]+Δtv[0]+Δt2M1f(x[1])

      The equation equiv to the following (where xM2=xTMx) => Optimization prob => Numerical schemes (Usually only conservative forces can use this energy)

      x[1]=argminF(x) for F(x)=12Δt2xx[0]Δtv[0]M2+E(x)

      Because: (applicable for every system; The first order der of F(x) reaches 0 for the min pt.)

      F(x[1])=1Δt2M(x[1]x[0]Δtv[0])f(x[1])=0x[1]x[0]Δtv[0]Δt2M1f(x[1])=0
      The Optimization Problem
      Newton-Raphson Method

      Solving optimization problem: x[1]=argminF(x) (F(x) is Lipschitz continuous)

      Given a current x(k) we approx the goal by 0=F(x)F(x(k))+F(x(k))(xx(k)) (Taylor Expansion)

      (For 2D: 0=F(x)F(x(k))+F2(x(k))x2(xx(k)))

      Steps:

      Newton’s Method finds extremum, but it can be min or max => finds 2nd order derivative (at local min: F(x)>0; at max: F(x)<0)

      image-20211130210036973

      For a function which second order derivative always larger than 0 (F(x)>0 everywhere) => F(x) has only one minimum

      Simulation by Newton’s Method

      For simulation: F(x)=12Δt2xx[0]Δtv[0]M2+E(x)

      Steps:

      Spring Hessian

      Hessian matrix is a second order derivative (sim to 2F/x2) => if p.d. so that the ONLY min point is found and has no maximum (sufficient but not neccesary condition)

      H(x)=e={i,j}[2Exi22Exixj2Exixj2Exj2]=e={i,j}[HeHeHeHe]

      The matrix in the last part is 3N x 3N, every vertex is a 3D vector. The first He is at (i,i), the He in the first line is at (i,j), …

      Positive Definite: Dep on every He (3 x 3): @ Lec.2, P48

      He=kxijxijTxij2+k(1Lxij)(IxijxijTxij2) ;xij=xixj

      where kxijxijTxij2 & (IxijxijTxij2) are already s.p.d. (Proof by multiplying by vT and v at both sides)

      But (1Lxij) can be negative when xij<Le (when compressed)

      Conclusion: when stretched He is s.p.d. and when compressed He may not be s.p.d. => A may not be s.p.d. either

      A=1Δt2M+H(x)=1Δt2Ms.p.d.+e={i,j}[HeHeHeHe]may not be s.p.d.

      (for smaller Δt => more p.d., actually sim to explicit integration with smaller time step => more stable)

      Positive Definiteness of Hessian

      When a spring is compressed, the spring Hessian may not be positive definite => Multiple minima

      image-20211130233603769

      Only in 2D/3D: In 1D: E(x)=12k(xL)2 and E(x)=k>0

      Enforcement of P.D.

      Some linear solver may require A must be p.d. in AΔx=b

      Solution:

      Linear Solvers

      Solving AΔx=b

      Jocobi Method

      vanilla Jocobi method (α=1) has a tight convergence req on A: Diagonal Dominant

      Other Solvers
      After-Class Reading

      Bending and Locking Issues

      The Bending Spring Issue

      A bending spring offers little resistance when cloth is nearly planar, since its length barely changes.

      image-20211207150345778

      A Dihedral Angle Model

      Every vertex of the 4 will be under a force as a function of θ: fi=f(θ)ui

      image-20211207150527413

      Conclusion:

      N1=(x1x3)×(x1x4) ; N2=(x2x4)×(x2x3) (Cross prod tells the dir. (for normal should be normalized) )

      E=x4x3

      Force:

      Explicit Derivative is difficult to complete / No energy

      A Quadratic Bending Model

      2 assumptions: planar case (OK for cloth); little stretching (mainly bending)

      image-20211207154104161

      Energy function: (vector * matrix * vectorT) (I is 3x3 identity)

      E(x)=12[x0x1x2x3]Q[x0x1x2x3];Q=3A0+A1qqTR12×12;q=[(cotθ1+cotθ3)I(cotθ0+cotθ2)I(cotθ0cotθ1)I(cotθ2cotθ3)I]R12×3E(x)=3qTx22(A0+A1); and E(x)=0 when triangles are flat (co-planar)

      Actually finding the laplacian (/ curl), when flat, no curl => E(x)=0

      Q is a constant matrix => The function is quadratic

      Pros & Cons

      The Locking Issue

      In the mass-spring / other bending models, assuming cloth planar deformation and cloth bending deformation are independent

      In zero bending case: LHS - OK; RHS - For stiff spring NO (Locking Issue) => short of DoFs.

      image-20211207161346501

      For a manifold mesh, Euler formula: #edges = 3#vertices - 3 - #boundary_edges

      If edges are all hard constraints: DoFs = 3 + #boundary_edges

      => no perfect solution

      Shape Matching

      ..

       

      Lecture 6 Constraint Approaches: PBD / PD / …

      Strain Limiting and Position Based Dynamics

      The Stiffness Issue

      Real-world fabric resist strongly to stretching

      But in simulation, only increasing the stiffness can cause: (More expensive / time-costing)

      A Single Spring

      If a spring is infinitely stiff -> Length = const

      Def a constraint func: ϕ(x)=xixjL=0 (Rest length L)

      image-20211207171930711

      Def a proj func: suppose in a R6 space for a pos of x want to move into a rational area (in blue, with boundary δΩ) with a shortest path: want the projection

      image-20211207172245723

      {xinew ,xjnew }=argmin12{mixinew xi2+mjxjnew xj2}such that ϕ(x)=0{xnewProjection(x)xinew ximjmi+mj(xixjL)xixjxixjxjnew xj+mimi+mj(xixjL)xixjxixjϕ(xnew )=xinew xjnew L=xixjxi+xj+LL=0

      The 2 new points’ substraction equals the original length => satisfying the constraint

      By default mi=mj, but can also set mi= for stationary nodes (can be just ignored)

      Multiple Springs

      Gauss-Seidel Approach

      Approach every spring sequentially in a certain order (needs a lot of iter to converge)

      image-20211207173051472

      image-20211207173152018

      Jocabi Approach

      Projects all edges simultaneously (good for parallization) and then linearly blend the results

      image-20211207173546431

      sum up all and update together with a weighted average.

      Position Based Dynamics (PBD)

      Based on the proj func, similar to shape matching

      image-20211207174356898

      Properties:

      Pros & Cons:

      (Usually in 3D softwares)

      Strain Limiting

      Some normal updates + strain limiting (corrections)

      image-20211207203039808

      e.g.: The constrain is not strictly equal to rest length L but some constraints (σ - stretching ratio as a limit)

      image-20211207203146606

      xnewProjection(x)σ1Lxixjσmin(max(σ,σmin),σmax)[σmin,σmax]{xinew ximjmi+mj(xixjσ0L)xixjxixjxjnew xjmjmi+mj(xixjσ0L)xixjxixjPBD: σ1;No limits: σmin,σmax

      Can be used to simulate some cloth whose stiffness increases when being stretched heavily.

      The constraints can be used for reducing oscillation for FEM / …

      Triangle Area Limit

      Limit the triangle area. Define a scaling factor first s=min(max(A,Amin),Amax)/A (Mass center no change)

      {xinew ,xinew ,xknew }=argmin12{mixinew xi2+mjxjnew xj2+mjxknew xk2}

      image-20211207204422565

      image-20211207204545942

      Properties:

      Projective Dynamics

      Main Idea

      Instead of using projections to update x directly as in PBD, projective dynamics uses projection to define a quadratic energy (identical to spring energy, but have 2 intermediate results xe,inew & xe,jnew)

      E(x)=e={i,j}k2(xixj)(xe,inewxe,jnew)2=e={i,j}k2xixjxixjxixjLexixjxixj2=e={i,j}k2(xixjLe)2{xe,inewxe,jnew}=Projectione(xi,xj) for every edge e

      Finding the corresponding forces: (assumption: (xe,inewxe,jnew)=const) (identical to spring force)

      fi=iE(x)=ke: ie(xixj)(xe,inewxe,jnew)=ke: ie(xixjLe)xixjxixj

      => althought the results are the same, the Hessians are different

      image-20211214000539781

      E(x)=k2(x0x1)(x0,0newx0,1new)2+k2(x0x2)(x1,0newx1,2new)2+k2(x0x3)(x2,0newx2,3new)2+k2(x1x2)(x3,1newx3,2new)2+k2(x2x3)(x4,2newx4,3new)2
      f0=[k(x0,0new x0,1new )k(x0x1)]+[k(x1,0new x1,2new )k(x0x2)]+[k(x2,0new x2,3new )k(x0x3)]
      H=[3kIkIkIkIkI2kIkIkIkI3kIkIkIkI2kI]

      The Hessian is the 2nd order derivative of energies to the vertices. The coefficients (of the diagonal) can be regarded as the numbers of the springs the vertices connected (3, 2, 3, 3 for vertex 0-3)

      This matrix is simple and constant A=1Δt2M+H

      => Applying implicit integration

      Simulation by PD

      Combination fo the projectiive dynamics with implicit time integration. Use a direct LU solver and can factorize A once (usually the most expensive part is the factorization and here the costs can reduce)

      image-20211214002319156

      Principles

      Newton’s Method calculate H as Hessian

      (1Δt2M+H)HessianΔx=1Δt2M(x(k)x[0]Δtv[0])+f(x(k))Negative gradient

      But in practice, H usually not calculated exactly (The performance dep. on how well the Hessian is approx., identity matrix can even be used (but convergence takes longer))

      Pros and Cons

      Pros

      Cons

      Constrained Dynamics

      This method is suitable for strong constraints / contacts / articulated rigid bodies (ragdoll animation) …

      A critical problem: what if constraints/forces are very stiff (or even infinitely stiff)

      example: multiple rigid bodies: human body can be regarded as several rigid bodies with very strong constraints (joints)

      Compliant constraint: ϕe(x)=xeixejLe

      The energy can be rewrited:

      E(x)=e12k(xeixejLe)ϕe2=12ΦT(x)C1Φ(x)f(x)=E=(EΦΦx)T=JTC1Φ=JTλ

      Let N be the number of vertices and E be the number of constraints (C - compliant matrix, J - Jacobian, λ - Lagrangian multipliers (Dual variables))

      Φ(x)RE ,[ϕ0ϕ1ϕE] ;C=[1k1k]RE×E ;J=ΦxRE×3N ;λ=C1ΦRE

      By implicit Integration: (2 sets of var: primal var x (or v=x˙) and the dual var λ)

      {MvnewΔtJTλnew=MvCλnew=ΦnewΦJ(xnewx)ΦΔt Jvnew[MΔtJTΔtJC][vnewλnew]=[MvΦ]

      Infinite stiffness in this case could be solved by C0 (smaller C makes the system better)

       

      Lecture 7 Linear Finite Element Method (FEM)

      Linear Finite Element Method (FEM)

      The Linear FEM Assumption

      In a nutshell, linear FEM assumes that for any point X in the reference triangles, its deformed correspondence is x=FX+c (“Linear” for this linear formula)

      image-20211214162153805

      For any vectors between 2 points we can use F to convert it from ref to deformed state:

      xba=xbxa=FXb+cFXac=FXba

      For example for edge 10 & 20: (indep on translation)

      FX10=x10 ; FX20=x20F=[x10x20][X10X20]1

      Problem: F is related to deformation but contains rotations (want to also reduce rotation)

      Green Strain

      Ideally, we need a tensor to describe shape deformation only. Recall that SVD gives F=UDVT, where only VT and D are relevant to deformation.

      image-20211109165409060

      Get rid of U as: G (Green Strain, symmetric, only 3 var) (FTF=VDUTUDVT)

      G=12(FTFI)=12(VD2VTI)=[εuuεuvεuvεvv]

      Strain Energy Density Function

      Consider the energy density per ref area: W(G):

      Total energy:

      E=W(G) dA=ArefW(εuu,εuv,εvv)

      The Saint Venant-Kirchhoff Model (StVK): (λ and μ are Lamé parameters)

      W(εuu,εvv,εuv)=12(εuu+εvv)2+μ(εuu2+εvv2+2εuv2)
      WG=[Wεuu12Wεuv12WεuvWεvv]=[2μεuu+λεuu2μεuv2μεuv2μεvv+λεvv]=2μG+λtrace(G)I=S

      (Trace is the diagonal sum; S - Second Piola-Kirchhoff stress tensor (some form of force density))

      Force:

      fi=(Exi)T=Aref(Wxi)T=Aref(Wεuuεuuxi+Wεvvεvvxi+Wεuvεuvxi)T

      The Wεij are solved before, we concentrate more on the next part

      Recall that: F=[x10x20][r10r20]1=[x10x20][abcd]=[ax10+cx20bx10+dx20]

      By def

      G=12(FTFI)=[12(ax10+cx20)T(ax10+cx20)1212(ax10+cx20)T(bx10+dx20)12(ax10+cx20)T(bx10+dx20)12(bx10+dx20)T(bx10+dx20)12]
      {εuux1=a(ax10+cx20)Tεvvx1=b(bx10+dx20)Tεuvx1=12a(bx10+dx20)T+12b(ax10+cx20)Tεuux2=c(ax10+cx20)Tεvvx2=d(bx10+dx20)Tεuvx2=12c(bx10+dx20)T+12d(ax10+cx20)T

      (the method is too complex => simplification)

      f1=Aref(Wεuuεuux1+Wεvvεvvx1+Wεuvεuvx1)T=Aref(Wεuua(ax10+cx20)T+Wεvvb(bx10+dx20)T+Wεuv12a(bx10+dx20)T+Wεuv12b(ax10+cx20)T)=Aref([ax10+cx20bx10+dx20][Wεuua+Wεuv12bWεuv12a+Wεvvb])=Aref[ax10+cx20bx10+dx20]deformation gradient[Wεuu12Wεuv12WεuvWεvv]Second Piola-Kirchhoff Stress[ab]=ArefFS[ab]
      [f1  f2]=Aref FS [r10  r20]1 ;f0=f1f2

      For tetrahedron (3D ref -> 3D deformation)

      Finite Volume Method (FVM)

      For simple triangle / tetrahedra, FVM is equivalent to linear FEM

      FVM considers force calculation in an integration perspective other than a differentiation perspective

      Traction and Stress in FVM

      Basic Ideas in FVM

      General Case:

      The Triangle Case: Take the midpoints which will be on the curve to estimate (the force contributes the same to the 3 vertices, the integral is for x0)

      image-20211221120134281

      In 3D (tetrahedra): FVM does the surfacial integral

      image-20211221170657979

      Different Stresses

      The stress tensor here: mapping from the interface normal to the traction (but kind of different from FEM)

      image-20211221171303566

      Output \ InputInterface normal N in the ref state (Undeformed)Interface normal n in the current state (Deformed)
      Traction in the ref state (unformed)Second Piola-Kirchhoff stress (S)-
      Traction in the cur state (formed)First Piola-Kirchhoff stress (P=FS)Cauchy Stress (σ=det1(F)PFT=det1(F)FSFT)

      (F - Deformation gradient)

      The Finite Volume Method

      image-20211221214405998

      f1=σ6(x01×x21+x21×x31+x31×x01)   (Replace σ with P (=> normal))=P6(X01×X21+X21×X31+X31×X01)(The term including X is constant => precomputed)=FS6b1(b1 is the constant)

      Second Piola-Kirchhoff Stress: S=WG as in previous FEM formulation

      => To find b1 (In the figure below, X1 will be computed via b1)

      image-20211221220647936

      The dot product of the edge 10 with X1

      X10Tb1=X10T(X01×X21+X21×X31+X31×X01)=X10T(X21×X31)=X01T(X31×X21)=6 Vol

      Meanwhile: X20T & X30T equal to 0

      [X10X20X30]Tb1=[6 Vol00] ;[X10X20X30]Tb2=[06 Vol0] ;[X10X20X30]Tb3=[006 Vol][b1b2b3]=6 Vol [X10X20X30]T=1det([X10X20X30]1)[X10X20X30]T

      Algorithm (Explicit)

      image-20211221225719287

      -> Various material (Heterogeneous materials) / various elastic / …

      Hyperelastic Models

      StVK: cannot handle flips / … (only widely used in graphics)

      Hyperelastic: More used in engineering (material / mechanics /… ) => more general

      Isotropic Materials

      Isotropic

      Force on any direction affects the same (with same behaviors) (-> ansotropic)

      F - Green Strain -> SVD; D=λ0,λ1,λ2 - Principal stretches: the singlar value of F

      P(F)=P(UDVT)=UP(λ0,λ1,λ2)VT

      In many lit, parameterize P(IC,IIC,IIIC) by principal invariants: (C=FFU is the right Cauchy-Green deformation tensor)

      IC=tr(C)=λ02+λ12+λ22 ;IIIC=tr(C2)=λ04+λ14+λ24IIC=12(tr2(C)tr(C2))=λ02λ12+λ02λ22+λ12λ22

      Models

      Solving for First PK Stress (P)

      Use the principal stretch for computation => first P-K stress

      P(λ0,λ1,λ2)=[Wλ0Wλ1Wλ2]P=UP(λ0,λ1,λ2)VT

      Algorithm

      Note: for the same model the first PK in this algorithm should equiv to P=FWG

      image-20211222113045756

      The Limitation of StVK

      The compression (with high stress) and flip cannot be handled with StVK model (even get stable in the flipped side) => UNSTABLE

      => Neo Hookean solve the problem

      image-20211222113516604

      Lecture 8 Linear Finite Element Method II

      Nonlinear Optimization

      Gradient Descent

      Solving x=argminF(x) => The negative gradient dir is the fastest descent dir

      image-20211222114523262

      Algorithm

      image-20211222114657200

      Step Size Adjustment

      Linear search methods: Exact line search (solve another optimization problem) & Backtracking line search (β is a smaller num (~ 0.3-0.4) -> Wolfe Condition)

      image-20211222114751055

      Descent Directions

      The dir d(x) is descending if a sufficiently small step size α exists for: F(x)>F(x+αd(x)) (at least descending)

      ie. F(x)d(x)>0 (the same side with the negative dir)

      image-20211222211524599

      With line search. use any search dir as long as it’s descending F(x(0))>F(x(1))>F(x(2))>F(x(3))>

      image-20211222211906481

      Descent Methods

      A unified descent framework:

      image-20211222212628550

      Differences between these methods:

      image-20211222212648030

      Total cost = Per-iter cost * Num of iter

      image-20211222212838495

       

      Lecture 9 Collision Handling

      Collision Detection

      Collision Detection Pipeline

      image-20211228113838360

      Board-Phase Collision Culling

      Spatial Hashing (Partitioning)

      Static space division

      Spatial partitioning divides the space by a grid and stores objects into grid cells

      image-20211228114157341

      Cell 12: t4Cell 13: t4Cell 14: t5Cell 15: t5
      Cell 8: t2Cell 9: t3Cell 10: t3, t5Cell 11: N/A
      Cell 4: t0, t2Cell 5: t0, t3Cell 6: t3Cell 7: t1
      Cell 0: t0Cell 1: N/ACell 2: t1Cell 3: t1

      For example for t3: Cell 5, 6, 9, 10 => other triangles stored in these cells: t0 & t5 => pairs: t0t3 & t3t5 and do further detection

      For moving objects, just expand the object region (add more)

      image-20211228114759191

      Optimization for memory costs: Memory wastes in empty cells and too many meshes should be built in 3D. Instead of allocating memories to cell, build an object-cell list and sort them

      e.g., (not relavent to space)

      Then sort by cell ID => easy to find out what is stored in each cell

      Morton Code

      One question is how to define the cell ID. Using the grid order is not optimal, since it cannot be easily extended and it is lack of locality (LHS). Morton code uses a Z-pattern instead (RHS).

      image-20211228164349900

      Bounding Volume Hierarchy (BVH)

      Bounding volume hierarchy is built on geometric/topological proximity of objects.

      External Object

      image-20211228170551965

      The parts that are not like to collide in different boxes => test if box intersection occurs

      To find elem potentially in collision with an obj, just traverse the tree:

      image-20211228171141111

      Internal Intersections (Self Collision)

      2 procedures in the tree

      BVH Models

      Axis-aligned bounding box (AABB) / Spheres / Oriented bounding box (OBB)

      Two AABBs intersect if and only if they intersect in every axis; shperes: need to compute the distance

      image-20211228172043462

      Narrow-Phase Collision Test

      Discrete Collision Detection (DCD)

      Actually not collision detection but intersection detection

      DCD tests if any intersection exists in each state at discrete time instant: x[0],x[1],

      To a triangle mesh, the basic test is edge-triangle intersection test: easy and robust

      image-20211228173830786

      (The t means the interpolation value between the points other than time)

      Tunneling

      The difference between collision and intersection detection:

      image-20211228174514285

      Problem: The green triangle at state x[0] and x[1] doesn’t have intersection with the blue one but actually the position relationship changes => tunneling problem (penetrating through each other)

      Not usually occurs for bulky objects but can happen for thin surfaces

      Continuous Collision Detection (CCD)

      CCD tests if any intersection exists between two states: x[0] and x[1]

      To a triangle mesh, there two basic tests: vertex-triangle and edge-edge tests

      Vertex-Triangle:

      image-20211228174836931

      (x30=x3x0 which is a function about time)

      Need to solve a third order equation (with one variable) => has formulae solution (not recommended, need to solve cubic root) / use binary division (0, 1) and get the first sol

      Edge-Edge:

      image-20211228212022858

      Issues with CCD:

      Continuous Collision Response Approaches

      The Two Continuous Collision Response Approaches

      Given the calculated next state x[1], we want to update it into x[1], the path from x[0] to x[1] is intersection-free

      image-20211229005436625

      Interior Point Methods (Cont.)

      Log-Barrier Interior Point Methods

      For simplicity, just consider the Log-barrier repulsion between 2 vertices (usually in graphics do a cutoff at some distance)

      image-20211229010338910

      E(x)=ρlogxijfi(x)=iE=ρxijxij2 ;fj(x)=jE=ρxijxij2

      Implementation

      Formulate the problem:

      image-20211229010911256

      x[1]argminx(12xx[1]2ρlogxij)

      Gradient Descent:

      image-20211229011306261

      The step size α must be adjusted to ensure that no collision happens on the way. To find α, we need collision tests.

      Impact Zone Optimization (Cont.)

      The goal of impact zone optimization is to optimize x[1] until it becomes intersection-free. (This potentially suffers from the tunneling issue, uncommon)

      x[1]argminx12xx[1]2

      Such that:

      {C(x)=(x3b0x0b1x1b2x1)N0 For each detected vertex-triangle pair C(x)=(b2x2+b3x3b0x0b1x1)N0 For each detected edge-edge pair 

      image-20211229012241218

      The optimizations are always trying to get close to the solution

      Rigid Impact Zones

      The rigid impact zone method simply freezes vertices in collision from moving in their pre-collision state. It’s simple and safe, but has noticeable artifacts

      Only used when all other methods not work

      A Practical System

      image-20211229012717208

      Untangling Cloth (Discrete)

      Intersection Elimination

      Consider how to eliminate existing intersections, but without using any collision history

      image-20211229013144634

      Useful when there are already intersections in simulation, due to:

      In this case, we don’t require the simulation is to always intersection-free

      Eliminating cloth-volume and volume-volume intersections is straightforward: simply pushing vertices/edges in the volume out

      image-20211229212934397

      Untangling Cloth

      Cloth-cloth intersection is complicated (don’t have a clear definition of inside and outside)

       

      Lecture 10 Surface Waves

      Two Types of Simulation Approaches

      image-20220104125629085

      A Height Field Model

      Height Field

      In 2D, a (1.5D) height field is a height function h(x) (cannot have 2 h at one position) and the velocity is also a function u(x) (negative -> opposite direction)

      image-20220104165850951

      The height field and velocity can be represented by: (hu means the volume flowing through a position at the unit time)

      {dh(x)dt+d(h(x)u(x))dx=0du(x)dt=u(x)du(x)dxadvection1ρdP(x)dx+a(x)ext.

      image-20220104170150684

      Ignoring advection and external acceleration, we get a simple form of velocity derivative: (ρ - density; P(x) - pressure)

      du(x)dt=1ρdP(x)dx

      From x to x+dx: dP(x)=P(x+dx)P(x) => the velocity changes due to the change of pressure (also dep on density)

      u(x) increases when dP(x)<0 and increases when dP(x)>0

      Shallow Wave Equation

      Simplification of the two equations by assuming shallow wave (the wave is small and the height change can be neglected)

      {dhdt+d(hu)dx=0dhdt+udhdx+hdudx=0dudt=1ρdPdxdtdx{d2hdt2+hd2dxdt=0d2udxdt=1ρd2Pdx2

      Eliminate u and formulate the shallow wave equation (indep on the velocity field)

      d2hdt2=hρd2Pdx2

      Finite Differencing

      Discretization

      image-20220104173532131

      Finite Difference

      Use the difference to approximate the derivative

      image-20220104173745771 image-20220104173759400 image-20220104174219237

      Second-Order Derivatives (central differencing in the fluid field)

      Applying central differencing twice to estimate height d2hi/dt2 (Laplacian)

      d2hi(t0)dt2dhi(t0+0.5Δt)dtdhi(t00.5Δt)dtΔthi(t0+Δt)+hi(t0Δt)2hi(t0)Δt2

      Also for pressure d2P/dx2

      d2Pidx2dPi+0.5dxdPi0.5dxΔxPi+1+Pi12PiΔx2

      Discretized Shallow Wave Equation

      Want to obtain the height field in the next time step

      d2hdt2=hρd2Pdx2hi(t0+Δt)=2hi(t0)hi(t0Δt)+Δt2hiΔx2ρ(Pi+1+Pi12Pi)

      Problem: Volume loss or generation

      Volume Preservation

      We want the volume to stay the same

      Suppose hi(t)=hi(tΔt)=V, but by using the formuation in the prev section, we will get

      hi(t+Δt)=2hi(t0)Vhi(t0t)V+Δt2hiΔx2ρ(Pi+1+Pi12Pi)=V+Δt2hiΔx2ρ(Pi+1+Pi12Pi)

      and we cannot guarantee Δt2hiΔx2ρ(Pi+1+Pi12Pi) is 0

      Solution 1

      Modify scheme into: (not using hi as the coefficient => using hi1+hi2 and hi+1+hi2 instead)

      hi(t0+Δt)=2hi(t0)hi(t0Δt)+Δt2hiΔx2ρ(Pi+1+Pi12Pi)hi(t0+Δt)=2hi(t0)hi(t0Δt)+Δt2Δx2ρ((hi1+hi2)(Pi1Pi)+(hi+1+hi2)(Pi+1Pi))

      Can be understanded as the water exchanges between hi and hi+1. The water flowing from the LHS will flow into the RHS => Volume preserved

      Solution 2

      Assuming hi in the right term is const

      hi(t0+Δt)=2hi(t0)hi(t0Δt)+Δt2HΔx2ρ(Pi+1+Pi12Pi)So that hi(t0+Δt)=V+Δt2HΔx2ρ((Pi+1Pi)+(Pi1Pi))Must be 0

      Pressure

      The pressure is related to the water height: Pi=ρghi

      hi(t0+Δt)=2hi(t0)hi(t0Δt)+Δt2HgΔx2α(hi+1(t0)+hi1(t0)2hi(t0))

      Replace the term with a constant α

      Viscosity

      Viscosity tends to slow down the waves (Add a viscosity constant β: too high: not viscos; too low: not like fluid)

      hi(t0+Δt)=hi(t0)+β(hi(t0)hi(t0Δt))Momentum here+α(hi+1(t0)+hi1(t0)2hi(t0))

      Boundary Conditions

      Algorithm with Neumann Boundary

      image-20220105172628214 – to 3D – image-20220105172707786

      Two-way Coupling

      The coupling will be two ways: between liquid->solid and solid->liquid (also can be bubbles / soft solids / …)

      Key: how to expel water out of the gray cell regions

      image-20220105173316343 image-20220105173333356

      Virtual Height

      Idea: Set up a virtual height vi so that hireal_new=hiei

      image-20220105173538947

      {Left:  hiei=hi+β(hihiold )+α(νi+1+hi+1+hi12νi2hi)=hinew +α(νi+12νi)Right:  hi+1ei+1=hi+1+β(hi+1hi+1old )+α(hi+2+νi+hi2νi+12hi+1)=hi+1new +α(νi2νi+1)

      Poisson’s Equation

      The outcome is Poisson’s equation, with νi and νi+1 being unknowns

      2νiνi+1=1α(hinewhi+ei)=biνi+2νi+1=1α(hi+1newhi+1+ei+1)=bi+1

      Regard as a Laplacian operator for every equation (νi1 & νi+2 are all 0, no virtual height)

      [2112][νiνi+1]=[bibi+1][112112111]=[νi1νiνi+1νi+2]=[0bibi+10]

      Algorithm

      Process: Set boundary condition -> Mask (which grids need solving) -> PCG_Solve() to retrun ν -> Update the height field (Usually need to multiply a factor, prevent too large waves when dragging the cube <- explicit’s unstability)

      image-20220105174829562

      Rigid Body Update

      Estimate the floating (buoyancy) force by the actual water expelled in every column

      image-20220105211306589

      fi=ρgΔx(hihinew)in 2Dfi,j=ρgΔA(hi,jhi,jnew)in 3D

      Also need to consider rotation => torque